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We consider the Generalized Gibbs Ensemble (GGE) in the context of global quantum quenches in XXZ 
Heisenberg spin chains. Embedding the GGE into the Quantum Transfer Matrix formalism we develop 
an iterative procedure to fix the Lagrange-multipliers and to calculate predictions for the long-time limit 
of short-range correlators. The main idea is to consider truncated GGE's with only a finite number of 
charges and to investigate the convergence of the numerical results as the truncation level is increased. 
As an example we consider a quantum quench situation where the system is initially prepared in the 
Neel state and then evolves with an XXZ Hamiltonian with anisotropy A > 1. We provide predictions 
for short range correlators and gather numerical evidence that the iterative procedure indeed converges. 
The results show that the system retains memory of the initial condition, and there are clear differences 
between the numerical values of the correlators as calculated from the purely thermal and the Generalized 
Gibbs ensembles. 
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1 Introduction 

Thermalization is a physical process whereby an isolated system reaches thermal equilibrium through 
the interactions of its constituents. Energy is conserved by time evolution, therefore the temperature of 
the equilibrated system is determined solely by the available energy at the beginning of the relaxation 
process. 

In quantum mechanics time evolution is represented by a unitary operation on the Hilbert space. 
When at t = the system is in a state ^>q, the density matrix evolves according to 

cn ; Pit) = |*(t)X*(*)| = e- im \^o){^o\e im 

■ * ■ Therefore, the density matrix itself never thermalizes. Instead, thermalization in quantum mechanics 

means that the expectation values of physical observables (or equivalently the matrix elements of the 
reduced density matrices) assume values which are equal to those calculated in a thermal ensemble. In 
a generic quantum mechanical system this is expected to happen, if there are no external driving forces 
present [TJ|2]. 

The situation is different in the case of integrable models, in which there exists an infinite family of 
mutually commuting quantum charges such that the Hamiltonian is a member of the infinite series: 

---' [Qj,Qk]=0 for j,fc = l...oo, Qi~H (1.1) 

It follows that time evolution conserves the expectation values of all charges, which prevents thermalization 
in the usual sense. Instead, relaxation to a Generalized Gibbs Ensemble (GGE) was proposed in [3J. The 
main idea is to construct a thermodynamic ensemble where the statistical weights depend also on all 
the higher charges and not only on the energy. It was proposed that the long-time limit of quantum 
observables after a quantum quench should be described by a GGE with appropriately chosen Lagrange 
multipliers for the individual charges. 

A huge amount of work was dedicated to study global quantum quenches in integrable models (see 
[H [5] and references therein) . Most of the available results concern theories equivalent to free fermions 
[3 [3 El El HOI HU HH US] or interacting models in the CFT limit [H]. One of the most important goals 



is to determine whether the GGE provides a valid description of the possible steady states arising in 
non-equilibrium situations. In the case of free theories the answer seems to be a general affirmative yes 
|15ll51 IT^[Tffl[TTl[T5l[TT?ll2Ull2I) . whereas the question is still unresolved for genuinely interacting theories. 
A further open problem in interacting theories is how to use the GGE to make actual predictions for 
physical observables. 

In Bethe Ansatz solvable models with only a single particle type the GGE hypothesis leads to a 
generalization of the Thermodynamic Bethe Ansatz (TBA) framework [2121 123] . Originally the TBA was 
developed to treat the purely thermal case (24] [25], but the addition of the higher charges does not modify 
its essential properties [23 . In principle the generalized TBA (or GTBA) encodes all thermodynamic 
properties and it can give predictions for the correlation functions as well [261 1271 I28j . In the Lieb-Liniger 
model studied in [22, 25jS], however, this can not be carried out because the construction of the higher 
charges seems to be an insurmountable problem |29| and therefore the Lagrange multipliers entering the 
GGE can not be fixed. Instead, information about the overlaps between pre-quench and post-quench 
states was used in [22] to set up the GGE. 

It is a very natural idea to apply the GTBA formalism to XXZ spin chains, where all higher charges 
can be constructed using their definition via the transfer matrix, or alternatively with the help of the 
boost operator [3U] EH] • However, the infinite family of particle types (the strings of different length) lead 
to an infinite set of GTBA equations. Although it has been demonstrated that even such an infinite set can 
be treated numerically [32], fixing the Lagrange multipliers for the GGE through TBA seems extremely 
difficult. 

An alternative to the TBA is the Quantum Transfer Matrix (QTM) method [321 GH], which leads to 
a single non-linear integral equation replacing the infinite set of the TBA [35| . Moreover, there are very 
efficient methods available to compute short-range correlations with the QTM [3(51 [371 [35] • The simplicity 
of the QTM with respect to the TBA makes it an excellent candidate to describe the GGE. The idea of 
adding higher charges to the QTM is not new: it was used in [3J3 [3D] to calculate thermal conductivities 
and in [411 I42j to study the phase diagram and thermal properties of integrable spin chains with competing 
interactions. 

In the present paper we apply the QTM formalism to set up the GGE for the spin chains. As an 
application we consider a specific example of a global quench, namely time evolution starting from the 
Neel state (in terms of the anisotropy A, this corresponds to a quench from A = oo to finite A). Although 
the relaxation of the antiferromagnetic order was already studied by different methods in [131 IS], these 
works did not consider the long-time limit of correlation functions. 

The paper is organized as follows. In Section 2 we pose the problem in general terms and explain our 
procedure to obtain the GGE as a limit of truncated GGE's. In Section 3 we present the extension of the 
QTM method to include a finite number of higher charges. Section 4 includes calculations concerning the 
explicit form of the conserved charges and their mean values in the Neel state. In Section 5 we present 
our numerical results and finally we conclude in Section 6. 

2 Global quenches in integrable models 

Consider a ID lattice model of L sites with periodic boundary conditions. The local spin variables are 
(Tj with j = 1 ... L and they can take K values, therefore the Hilbert-space of the system is 

U = ® L (C A ') . 

Consider a family of Hamiltonians of the form 

L 

3=1 

where it is understood that Uj and hjj+i are given by a translation of one-site and two site-operators u\ 
and h\ : 2 which might depend on a finite set of coupling constants. The generalization to Hamiltonians with 
multi-site interactions is straightforward. Note that in (|2.f p periodic boundary conditions are assumed. 

In this work we consider the situation of a sudden global quench. At t = we prepare the system 
in a state \^f(t = 0)) = l^o) which might be the ground state of a Hamiltonian Hq or a state prepared 
according to a well-defined rule. In the latter case we only require that \^o) be defined in a natural way 
for any L or at least any even L. 



At t = we change the parameters of the system such that time evolution for t > will be governed 
by the post-quench Hamiltonian H: 



|*(i)> = e- iHt \*c 



Consider a localized quantum observable O. In the examples to be investigated below O will be a 
correlation function of local spin operators on neighboring sites or only a few sites apart. The time 
dependence of this observable is 

0(L,t) = (*(t)|0|*(t)) - (*o|e 4jFft Oe-' im |* ). 

In the formulas above it was understood that all quantities depend on the volume L. We will be interested 
in the thermodynamic limit: 

0{t) = lim 0(L,t). 

L— >oo 

Moreover we consider the infinite time limit of the observable: 

6 = lim 0(t). (2.2) 

t— >oo 

In a generic non-integrable system we expect the phenomenon of thermalization. This means that all 
long-time averages are described by a thermal ensemble with an effective temperature T = 1//3: 

6 = ^Tr (Oe~ pH ) Z = Tr e~ &H . (2.3) 

Z 

Time evolution conserves the expectation value for the energy, therefore the Lagrange multiplier /3 can 
be fixed (at least in principle) by the requirement 

(H) = --^logZ=(* \H\* ). (2.4) 

A central statement of thermalization is that there is a single effective temperature T for all quantum 
observables. Note that in formulas (|2.3[) - (|2.4p an implicit infinite volume limit is understood. 

The situation is expected to be different when H is integrable, in which case the thermalization 
hypothesis does not hold. The reason for this is the following. In integrable models there exists an infinite 
set of conserved charges Qj. j — 1 . . . oo which commute among themselves: 

[Qj,Qk]=0. (2-5) 

The Hamiltonian is a member of the series, typically H ~ Q2 (the first charge Q\ is usually the momentum 
operator). It follows from (|2.5p that the expectation values of all the charges are conserved in time. The 
Qj are constructed as sums of localized operators, therefore (12. 3p should apply if the system thermalizes. 
However, the thermal ensemble would typically yield mean values for the conserved charges which differ 
from those measured in the initial state, therefore (|2.3|) can not be valid and thermalization can not occur. 
To describe the long-time average of observables in integrable models the Generalized Gibbs Ensemble 
(GGE) was proposed in [3J. Setting aside convergence properties for a moment the hypothesis of GGE 
can be formulated as follows: There exists a set of couplings {j3j} such that the long-time stationary state 
of the system is described by the density matrix 

Pgge = ^^ exp I - ]T fcQj Z GGE - Tr exp I - ]T ftQj . (2.6) 

Physical quantities in this ensemble are given by 

O gge = Tr {Op GGE ) ■ (2.7) 

Expectation values of the charges are conserved in time, therefore the Lagrange multipliers can be fixed 
(at least in principle) by requiring 

(Q 3 )=-^rlogZ GGE = (y \Q 3 \y ), j = l...oo. (2.8) 

apj 



The GGE hypothesis states that for any localized operator the expectation value obtained by the equations 
(|2.6p - (|2.8p is equal to the long-time average: 

O gge = 6. (2.9) 

Again, a crucial statement is that there is a single set of /3j which determines all physical observables. 
Note that all conserved charges need to be added to (|2.6[) . otherwise the GGE can not be complete. 

Equations (|2.6p ~ (|2.8p completely specify the GGE. However, in practice it is impossible to fix all the 
Pj and the convergence of the infinite sum is also not guaranteed. In this work we propose to obtain the 
GGE as a limit of an iterative procedure, where at step k we only consider the first k charges. This idea 
also appeared in the very recent paper |13| which concerned the quantum Ising chain. 

To be specific, we define sets of parameters {(3 ■ , j = 1 . . . k} which generate the density matrices of 
truncated GGE's: 

p« = -Lexp ( £/f Q; I 3« =Trcxp ( - £/f Q; I . (2.10) 



zw 



i=i / \ j=i 



We call the number k the truncation level. The j3, are chosen such that 
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in logZ( fe ) = (* |Q J |* ), j = l...k. (2.11) 

For the quantum observables we obtain a series 

{k) =Tr(Op {k A. (2.12) 

The actual GGE average is then given by the limit 

O gge = Um O ik) . (2.13) 

This procedure provides a well-defined recipe to obtain the GGE averages, but we are faced with 
the following questions: Does the limit (|2.13p exist? Does it exist for all localized physical observables or 
maybe even non-local quantities like correlation lengths? How do the convergence properties depend on 
the Hamiltonian and the initial state? 

While we can not answer these questions in their full generality, we will present one non-trivial example 
in XXZ spin chains where numerical evidence shows that the limit (|2.13p exists for short-range correla- 
tion functions. This way the GGE indeed gives predictions which could be compared to experiments or 
independent numerical calculations. 

3 Thermal and Generalized Gibbs ensembles for the XXZ spin 
chain 

In this work we consider quantum quenches in spin- 1/2 XXZ chains. The Hamiltonian is given by 

L M 

H xxz (J,A,h) = JX(oX+i +*Kh+ A (°M+i - !)) + h z2°r ^ l ) 

3=1 3=1 

In the following we will only consider the zero-field case h — 0. Moreover we set J = 1 and restrict 
ourselves to the regime A > 0. 

In the present section we set up the general framework of the GGE, independent of the details of the 
initial state. Here we just assume that the mean values of the conserved charges can be calculated with 
exact methods. Details of how to do this will be provided in the next section, and we give a few comments 
about more general situations in the Conclusions. 

The XXZ model is integrable for arbitrary A. Its spectrum is given by the Bethe Ansatz [441 H51 H51 Wf\ 
and the higher conserved charges can be constructed using the transfer matrix formalism. 



Consider the It-matrix acting on C <g> C given as 



R(u) ' 



sinh(u + 77) 



/sinh(w + 77) \ 

sinh(u) sinh(?y) 

sinh(Ty) sinh(w) 
V sinh(u + rj)J 



(3.2) 



Here u is the spectral parameter and 77 is a complex number related to the anisotropy by A = cosh 77. 
The monodromy matrix is constructed as 

T(u) = L M (u)...L 1 (u), 

where Lj (u) are local Lax-operators given by 

Lj(u) = Rqj(u), 

and the index stands for the auxiliary spin space. In this space T{u) can be written as 

'A{u) B(uY 



T(w) - \C(u) D(u)J ' 

where A(u),B(u), C(u), D(u) are operators acting on the spin chain. 
The transfer matrix is given by the trace 

t(u) = Tr T(V) = A(u) + D(u). 

The i?-matrix satisfies the Yang-Baxter equation [151 [41?] which leads to the commutativity of the 
transfer matrices: 

[t(u),t(«)]=0. 

This property is used to define the conserved charges of the model. It can be shown that 

r(0) = U, 

with U being the translation operator on the chain and it can be considered as the first conserved charge: 
U — e 1 ® 1 , where Q\ is the momentum operator. The other charges are defined as logarithmic derivatives 
of the transfer matrix at u = 0: 

du 



Qj=( — ) l0 §^)- (3-3) 



It was shown in [50] that the Qj defined this way are local in the sense that they are given as sums of 
products of spin variables such that they only span a finite segment of the chain of length j. We note that 
using the normalizations above the second conserved charge is 

Q2 = . r H X xz- 

2 sinn 77 

Eigenstates of the commuting family of transfer matrices can be constructed using the Algebraic Bethe 
Ansatz [51j . We choose the reference state \F) to be the ferromagnetic state with all spins up. Bethe states 
are then constructed by acting with the i3-operators: 

M 
\n 1 ,..., m ) = l[B( H )\F). (3.4) 

Here the complex variables fij are the rapidities of the interacting spin waves. Such a state is an eigenstate 
of the family of commuting transfer matrices if the Bethe equations are satisfied: 

d(M II ■ w 7 = 1, (3-5) 

AA smh(/i i -fik-V) 

where 

d(u) = (F\D(u)\F) = ( ,7 h( " } )\ (3.6) 

\smh(7j + f]J J 



It can be shown that the states p. 41) are identical to the states constructed by the coordinate Bethc 
Ansatz. 

In the regime A > 1 we have 77 £ K and the usual one-string solutions of the Bethe equations are of 
the form /x = iX — r)/2 with Ael. On the other hand, when A < 1 we have r\ = 17 with 7 £ R and the 
usual one-string solutions are /i = A — i"f/2 with A£R. The isotropic case (A = 1) can be obtained by a 
limiting procedure, but we do not consider this here. 

In the following we will need the eigenvalues of the conserved charges on the Bethe states. The Algebraic 
Bethe Ansatz yields the following transfer matrix eigenvalue: 



t(u, Mm) 



M 

n 

3=1 



sinh(u — /ij — rj) 
sinh(w — jij) 



M 



•d(«)IJ 
3=1 



Taking the logarithmic derivative gives 

Qilfii,... 



,PN, 



= iJ2<ij(^k)j 



Mi, 



sinh(u — jij + r]) 
sinh(u — /ij) 



,/Ujv) 



with 



Qj(v) 




sinh(u 
log 



/i - rj) 



sinh(u — n) 



Here we assumed that L > j such that the terms originating from d(u) are all zero. 
For the second charge (proportional to the one-string energy) we obtain 



92 (m) = 



cosh(ju) cosh(/i + -q) 



sinh(ju) sinh(/i + rj) 
The formulas for the higher charges can be written in the form 

?jO) =Gj(xo)-Gj(x+), 
where 



x 



cosh(/i) 
sinh(/i) 



x + 



cosh(fi + rj) 
sinh(/i + rj) 



and the Gj (x) are polynomials satisfying the recursion 



(3.7) 



(3.8) 



G j+1 (x) = (x"-l)—G j (x), 



G 2 (x) 



(3.9) 



Thermal ensembles for the XXZ spin chains can be constructed using the so-called Quantum Transfer 
Matrix (QTM) formalism |33l IM] . The goal is to construct the density matrix 



1 



P 



exp(-^Q 2 ) Z = Trexp(-/3Q 2 ). 



The central idea is to write down the Trotter-Suzuki decomposition of exp(— fiQ-i): 

N 



exp(-/3Q 2 ) 



-> 



(r- 1 (0)r(-/3/A)) 



N 



(3.10) 



Here N is called the Trotter number and the two relations above become equalities in the N — > 00 limit. 
At finite N the trace 



Z n , l = Tt(t- 1 (0)t{-P/N)) 



N 



can be interpreted as a partition function of a six- vertex model with L vertical and N horizontal lines. 
The vertical lines correspond to the original spin spaces and the horizontal ones are auxiliary spaces. The 
partition function can be obtained alternatively by "quantizing" the system in the horizontal direction 
introducing the "Quantum Transfer Matrix" which acts on the auxiliary spaces. The L->oo limit of the 
partition function can then be obtained by the largest eigenvalue of the QTM. It is known that the QTM 
is gapped in the sense that the second largest eigenvalue is separated from the largest by a finite distance 



even in the N — > oo limit. This means that all thermodynamic properties will be determined by the single 
leading eigenstate of the QTM. This state can be found by the Algebraic Bethe Ansatz and its transfer 
matrix eigenvalues can be computed both at finite N and in the Trotter limit. For the details of this 
procedure we refer the reader to the review |34j . 

In the Trotter limit the partition function can be expressed as 

logZ = -/£ + ..., 

where the dots denote exponentially small corrections in L and the free energy density is given by 

f = f du sinh7ylog(l + a(q;)) 
J c 2ni sinha>sinh(u; + 77) 

Here a(A) is an auxiliary function defined on the complex plane which satisfies the nonlinear integral 
equation (NLIE) 

„ ,,. f duj sinh2ri \oe( 1 + a(uj)) 
loga(X)=-0q 2 (X)- - r-rr-, . , , x r. (3.12) 

B v ; ,HK ' J c 2ni sinh(A - uj + rj) sinh(A - u - 77) v ' 

The contour C in the equations above depends on A. For the A > 1 regime considered in the present 
work it can be chosen as a union of two straight line segments: 

C = [-in/2 + a, in/2 - a] U [in/2 - a, -in/2 + a], 

where a < r\/2 is an arbitrary parameteqj. Note that the first line segment runs upwards and the second 
runs downwards. 

Thermodynamic properties of the spin chain can be calculated by taking derivatives of the free energy 
with respect to the physical parameters. Correlation functions are also accessible to this method. In |52] 
multiple integral formulas were calculated for the localized correlation functions. These were later found 
to factorize, ie. they can be expressed as sums of products of simple integrals [53] so that numerical results 
can be produced in a very efficient way [3S1 [371 155] . 

It is a very natural idea to use the Trotter-Suzuki decomposition to construct the truncated GGE 
density matrices (|2.10p . In fact this method was already worked out in [39 where it was used to obtain 
thermal conductivities. 

The main idea is that for any finite k the decomposition 

k ( r fe _ bo- \ N 

exp(-Y,PjQj) * I 1 - ^ J -y J ^ +o(l/AT)j (3.13) 

holds. The right hand side can be obtained as 

1 - Ej ^f jQj = (r- l (0)) k ~ 1 t{ Ui )t{u 2 ) . . . r(«*_x) + o(l/JV), (3.14) 

with some appropriately chosen numbers Uj . For example in the case of k — 3 and non-zero 83 a solution 
is 

Ul V N 2N U2 V N 2N' 

Assuming that the numbers Uj are found the partition function can be expressed as 



^N.L,k 



Tr 



, fc-i 



(r- x (0))* T(uiM«j)...T(ttt-i) . (3.15) 



JY 



This is equivalent to a partition function of a six-vertex model, which can be quantized in the horizontal 
direction, leading to a modification of the original thermal QTM with a different set of inhomogeneities. 
We can assume that the analyticity properties necessary for the construction of the Trotter limit hold also 
in the modified problem, at least in a small neighborhood of the purely thermal case. Then the Trotter 
limit can be taken and it leads to 

\ogZ Lik = -f^L + ... (3.16) 



1 When a magnetic field is added, a has to be large enough such that the contour encircles all Bethe roots of the QTM. 



with 



,( fc ) = _ f du sinh7ylog(l + g( fc )(cj)) 
J c 2iri sinh u sinh(oj + rj) 



Here a*-' -' is the auxiliary function solving the modified NLIE 

m„> \—* „ ,,x /" dw sinh 2n loe(l + a^ k '(u)) 

logo (fc) (A)=- > &«*(A)- / - . , - .. ' B ; , , ,. v yy -. (3.18) 

to v y Z->^ J ^ / j c 2m s i n h(A - w 4- rj) sinh(A - w - ry) v y 

Note that the structure of the source term reflects the form of the truncated GGE density matrix (|3.13p . 
Expectation values of conserved charges in the truncated GGE can be obtained by 

m=L^. (3.19) 

Instead of taking the derivative of the formula (|3.17p it is useful to express /W with the functions 
a( k ) = l/a( fe ). It can be shown that they satisfy 

loga( fe )(A)=-V^ 7 (A) + [ £L , ^h2,l g(l + aW(c)) 
b v ; /L^^^o ^ > J c 2 m S mh(A - w + rj) smh(A - u - rj) v ; 



with 
where 



x 



qr(n)=G j (x )-G j (x-) ) (3.21) 

cosh(/i) cosh(/i — 77) 

sinh(/i) sinh(/i — 77) 

The polynomials Gj are defined by (|3.9p . The free energy can be expressed with 0^ as 

(fc)= / d^ sinh??log(l + aW(cj)) 
J c 2iri sinh w sinh(w — 77) 

Taking the derivative of (|3.20p and (J3.22I) with respect to /3j we introduce the functions 

(fc) m= 1 da^(A) 1 aqW(A) 

° J ' l j o( fe )(A) dpj a«(A) 3& ' 

They satisfy the linear equations 

a (fe) (A) = - a' (A) + f — sinh2?7 a ^ H (3.23) 

> { ' J Jc 2 ™ sinh (A - u + 77) sinh(A - w - ry) 1 + a« (w) ' 

Finally the conserved charges in the truncated GGE are given by 

(Q.\-L f — Sinh?/ ^ )(CJ) r3 241 

W^ lv i / c 27risinhu;sinh(w-r ? )l + a( fe )(a;)' K ' 

With this we have finished the construction of the truncated GGE using the QTM formalism. The 
remaining task is to calculate correlation functions in the truncated GGE. 

In the derivation of the multiple integrals of [S3] for the thermal correlations it is not necessary to 
know the exact position of the Bethe roots of the QTM; the only required information is that they can 
be surrounded by appropriate contours and that there exists an auxiliary function o(A) which encodes 
the positions of the Bethe roots through the equation o(A) = — 1 and which has a well-behaving Trotter 
limit. These conditions also hold for the truncated GGE's, at least in the neighborhood of a thermal case 
with a finite fa. In the multiple integrals of [S5] the auxiliary function (or the combination 1/(1 + a(A)) 
plays the role of a weight function, therefore the representation of (52J is valid also in the truncated GGE, 
provided that zeroes of (1 + o(A)) do not cross the contours. Based on continuity and symmetry arguments 
we expect that the Bethe roots will still be situated on the imaginary axis (when A > 1) or the real axis 
(when A < 1). Our numerical findings support this expectation in all cases we encountered (see Section 



5). We conclude that the multiple integrals of [52J are also valid in the truncated GGE's, at least in the 
cases considered in the present work. 

Factorized formulas for the correlators were developed in [53J. This factorization of the multiple in- 
tegrals depends only on certain algebraic properties underlying the construction of correlation functions 
and not on the particular physical parameters (finite temperature or finite size) of the problem at hand 
[541 I55] . It follows that all the formulas already available for the thermal case are also valid in the trun- 
cated GGE, provided the calculations are carried out using the auxiliary functions a^ (A) which satisfy 
(|3.18p instead of Q3.12p . We refrain here from replicating the necessary formulas and refer the reader to 
the original papers [531 I3H| ■ 

The method we described can be applied to construct the GGE within the QTM formalism and 
to calculate correlation functions in the GGE. The only input from the initial states is through the 
expectation values 

(Q 3 ) = (*o|Qil*o) 

which are used to fix the Lagrange multipliers. These mean values should be obtained exactly with use 
of explicit formulas for the Qj or by other methods. In the next section we show one example where this 
task can be performed relatively easily. 

4 Quench dynamics from the Neel state 

We consider a quench with the Neel state as initial state: 

|*o> = pV) = |+- + - + -...}■ 

This state is not translationally invariant, the correlation functions involving an odd number of spin 
operators (such as the magnetization (erf) itself) have an obvious position dependence. However, we 
expect that in the long time limit translational invariance of all correlation functions will be restored |43| , 
in particular 

lim (a*) = 0. 

We note that even though the magnetization relaxes to zero, traces of the original antiferromagnetic order 
are expected to show up in the long-time limit of two-point correlators a^a^ +1 and <x? cr| ', x . Our numerical 
result show that this is indeed the case. 

As an alternative to the Neel state we could also choose the initial states 

\* ) ± = ^=(\N)±\AN)), 

where \AN) = | — I 1 h . . . ) is the anti-Neel state. These states lead to correlation functions which 

are translationally invariant at all times. Moreover, the expectation values of the conserved charges are 
equal to the respective values in the Neel state, because in the infinite volume case considered here all cross- 
terms of the form (N\Qj\AN) vanish. We note that the cross terms do not influence the time-dependent 
correlations either, because 

(N\e iHt Oe- iHt \AN) =0 (4.1) 

at any finite t, given that O is a localized operator and the infinite volume limit is already performeqj. 
Therefore the predictions of the same GGE hold in all three cases, but it is enough to calculate the 
conserved charges in \N). 

The Neel state is not an eigenstate of the Hamiltonian, therefore the transfer matrix itself can not be 
used to calculate the mean values directly and the explicit form of the Qj is required. 

4.1 Conserved charges in the Neel state 

Higher conserved charges of XXX and XYZ spin chains were considered in the papers |30l I3"T] . The 
basic tool of these papers is the boost operator [5_SJ [SB] given (in the XXZ case) by the formal expression 



2 In finite volume there are non-zero contributions to (|4.1|) when the t is large enough so that quasiparticles can travel around 
the volume and therefore a complete shift of the two states is possible. 



It can be shown that the boost operator generates time derivative of the transfer matrix: 

[B,T(u)]~f(u). (4.2) 

It follows that conserved charges are generated recursively as 

[B,Qj] ~ Qj+i- 

The equation above is to be understood up to additive constants. 

The authors of [201 EH] considered the action of the commutator (|4.3p and derived a recursive procedure 
to obtain all terms of the charges. They were able to analytically solve the recursion in the XXX case, 
whereas in the general XYZ case (including the XXZ chains) they derived the explicit form of the charges 
up to Q e . 

In order to translate the results of [201 121] into our conventions we define operators Qf M as 

l[B,Qf M ] = Qf™ (4.3) 

with the first member being 

oo 
nGM \ "" i _x „x i „V „V i \ „z z \ 

j=—oo 

Up to overall phase factors these operators coincide with the ones given in [301 131] and they differ from 
the Qj used in the present work in additive and multiplicative normalization. The additive constants can 
be fixed by requiring that all mean values vanish in the reference state \F); this follows from the fact that 
in our normalization 

(F\t(u)\F) = 1. 

The multiplicative normalization can be deduced by working out the proof of (14.21) presented in [21] using 
our normalizations. We find the relation 

Qf M -(F\Qf™\F) 
^ j ~ 2(sinh^-i • [ ' 

An important basic result of [31] for the XYZ chain is that each conserved charge Qj is a sum of terms 
of the form 

[(...(o- n x a i2 ) x a i3 )--- x a lt _ t ] ■ a H . (4.5) 

Here (i\, ?2, • • • , ii) is a sequence of sites in increasing order such that I < j and ii — i\ < j. The notation 
a and a refers to vectors constructed out of rescaled Pauli-matrices. In the XXZ case they are 

&=(a x ,a v ,VAc7 z ) a=(VAa x ,VAa y ,a z ). 

In the Neel state only those terms are non-vanishing which only include a z operators. From all possible 
terms in (|4.5p we find such products only when 1 = 2, because any cross-product would necessarily involve 
at least one a x or a v operator. Therefore it is enough to collect the terms of the form afaf +l . These terms 
only appear in the even charges Q^j , which is consistent with the general statement that the mean values 
of the odd charges vanish in any parity invariant state. 

The paper [31 provides formulas up to Qg. In order to obtain explicit representations for the higher 
charges we implemented the iteration procedure (|4.3[) with the symbolic manipulation program form [57j . 
The formal commutation relation (|4.3p is only valid on an infinite spin chain, but it can still be used to 
calculate the charges. Defining the finite size operators 

3=1 

l-l 

Q l 2 = (of of + afaf + Aafaf) + 5>^ +1 + a^ +1 + Aa*a* +1 ) 

3=1 
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we implemented the recursion 

){B l ,Q< j }=Q l j+1 . (4.6) 



At each step of the iteration additional boundary terms arise which start to propagate towards the middle 
of the chain. However, if I is large enough then the "bulk" of Q ■ is not affected and the coefficients of the 
different terms of Qf M can be read off from the middle of the chain. We used this method to obtain all 
charges up to 0i2- For the sake of brevity we only present the terms relevant for the Neel state, and only 
up to Q s : 

3 

3 

Qf M = ^(56A + 16A 3 )ctK+i ~ ( 64A + 8A3 K a 3+2 + 24AfJ K+3 + ■ ■ ■ (47) 

3 

Qf M = ^(1504A + 1312A 3 + 64A 5 )(t j z ct| +1 - (2912A + 1376A 3 + 32A 5 )ct|^+ 2 + 



(2400A + 480A 3 )o-j : crj : + 3 - 720Aa| a 



z z 

3 a 3+± 



(4.8) 



The dots represent terms with at least two a x or a y operators. 
The Neel state mean values read (up to Q12) 

Q$ = -2A 

Of = -8A 

Q% = -(160A + 32A 3 ) 

Qg = -(7808A + 3584A 3 + 128A 5 ) 

Q10 = -(709120A + 517632A 3 + 62976A 5 + 512A 7 ) 

Q12 = -(103467008A + 103763968A 3 + 23973888A 5 + 1036288A 7 + 2048A 9 ), 

where we defined Qf = ((N\Qf M \N) - (F\Qf M \F))/L. 

For the sake of completeness we also present the ferromagnetic expectation values of Qj . Defining 

Qf = (F\Qf M \F)/L we obtain 

Q% = A 

Ql -2A 

Q- = 16A + 8A 3 

Qf = 272A + 416A 3 + 32A 5 

0f = 7936A + 24576A 3 + 7680A 5 + 128A 7 

0f 2 = 353792A + 1841152A 3 + 1304832A 5 + 128512A 7 + 512A 9 . 

We stress that eqs. (|4.9p represent the constant terms entering the r.h.s. of (|4.4[) . The ferromagnetic 
eigenvalues of the higher charges all vanish in our normalization. 

5 Numerical implementation 

For our numerical calculations we considered the regime A > 1 . In this regime the XXZ chain is gapped: 
there are two ground states which become degenerate in the L — > 00 limit, but there exists a finite gap 
between them and the next excited state. In the limit of A — > 00 the Hamiltonian turns into the classical 
Ising model, with its two ground states given by the Neel and anti-Neel states. Therefore quenching to 
a large but finite A can be considered a "small quench": we expect that all physical quantities predicted 
by the GGE will be close to their respective values in the Neel state. Departure from these mean values 
should show up for smaller A. For our actual numerical calculations we chose the values A = 2, 3, 4, 5. 
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We implemented the NLIE (|3.18[) and solved it numerically by simple iteration, which converged even 
when the higher charges were added. It is known that the iteration method becomes ineffective in the 
low-temperature regime, where /^A is large. In this case a different formulation of the NLIE can be set 
up [34 1; however this was not needed for our purposes. 

We also implemented the factorized formulas of [53] for the short-range correlators 

{(T x a 2 ) , {a 1 cr 3 } 

for a — z,x. We checked our programs by calculating these correlators in the purely thermal case and 
comparing them to the values published in [37] and we found complete agreement. 

We tested our numerics for the higher conserved charges by evaluating them in the purely thermal 
case and at low temperatures (j3j = for j > 3 and j3 2 — > oo). In this limit the thermal averages quickly 
approach the eigenvalues in the two ground states, because the Hamiltonian is gapped. On the other hand, 
the ground state mean values can be found by standard Bethe Ansatz calculations [2S] : 

00 p—vl^l 

lim (Q 2j )=- V (2n) 2 < J - 1) -. , . (5.1) 

/3 2 ->oo X J/ nfl'oo cosh(r]n) 

We compared our numerical results with the formula above and found convincing agreement. Note that 
the A — > oo limit of (15. ip reproduces the leading terms of (|4.8p . given that the normalization (|4.4p is 
used. This is a consequence of the fact that the states (| N) ± \AN))/y/2 become the two ground states in 
the A — > oo limit. 

As a further numerical verification we also evaluated the high-temperature limit of the mean values. 
At high-temperatures any product of spin operators tends to zero, because the Pauli matrices are traceless 
operators. Therefore only the constant terms of Q 2 j contribute and from (I4.4p we obtain 



(F\Q^\F) 
pT"o^ 2j/ 2(sinh?7) 2 i- 1 



lim (Q 2 



We compared our numerical results to those calculated from (|4.9p and found complete agreement^. 

In order to find the numerical values of the Lagrange multipl 
Newton- Raphson methocQ The matrix elements of the Jacobian 



(k) 

In order to find the numerical values of the Lagrange multipliers fi- we used the multi-dimensional 



d (Q m ) 

d/3 n 

are easily obtained by taking a further derivative of (13.2411 . The Newton- Raphson method converged if 
started from a purely thermal state with an appropriate j3 2 - 

We observed that the numerics became sensitive to the choice of the integration contour for the NLIE 
as we add the higher charges. There are always two numbers < a± and 0,2 < r]/2 such that every 
a chosen from the interval [ai,«2] gives the same answer for all quantities up to a desired numerical 
accuracy. However, the available interval for a shrinks with the growing number of the charges. This can 
be explained by the fact that the higher source terms for the NLIE are more and more singular, which in 
turn requires that the integration contour is far enough from the origin in order for the calculations to be 
stable. As a rule of thumb we used the values a — 0.95 x 77/2. 

In order to gather information about the Bethe roots of the QTM we calculated the function L(X) = 
(l + o(A)) along the imaginary axis using (|3.18p . Zeroes of L(X) determine the positions of the roots, which 
are known to be purely imaginary in the thermal case (at zero magnetic field) . In all our examples it was 
found that although the position of the roots changes considerably as compared to the thermal case, they 
always stay on the imaginary axis. This provides a strong justification for the validity of the NLIE (|3.18p . 

Having convinced ourselves that our numerics is stable, we computed the predictions for the short- 
range correlators. Our calculations were carried out with a total of n = 400 discretization points. When 
a was chosen from the stability interval, then the numerical values did not change up to 6 digits when 
we increased n to 600, or when we slightly changed a. We conclude that the results presented below are 
accurate up to the last digit. 

Tables [1] and [5] include our numerical results. The correlation functions are also plotted as a function 
of the truncation level in figures [T] and [21 



3 The high-temperature limit of (Q23} could be calculated also from (|3.23[) - (|3.24|) using lim^^o tt(A) = 1. 
4 This was suggested by Gabor Takacs. 
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Regarding the Lagrange multipliers our data is not sufficient to determine whether they are convergent 
as a function of the truncation level. Note that convergence of the /3j themselves is not required by 
physical principles, because they can not be measured. In fact higher charges can modify the lower j3j 
considerably, because they involve terms which are already present in the lower charges. We also calculated 
the cumulated coefficients 71 and 72 defined as 



k 

E 
1=2 



/r^ = -7i£^j+i 



72 V^ < T : <>' 



3 



J U J+2 



Our results at A 



5 are 


below: 












fc 


2 


4 


6 


8 


10 


12 


7i 


19.462 


18.632 


20.264 


20.108 


20.793 


20.774 


72 





0.23873 


0.22747 


0.28444 


0.28171 


0.30213 



These data are not sufficient to determine whether 71 and 72 are convergent. 



Regarding the correlators we observe the following: 

• For all A and all distances the z — z correlators increase, the x — x correlators decrease as a function 
of the truncation level. The physical reason for this is evident: The Neel state has maximal z — z 
(anti-)correlations and identically zero x—x correlations. Addition of higher conserved charges should 
reproduce more details of the initial state, leading to the observed behavior. 

• All correlators seem to converge for A = 4 and A = 5. In order to convince ourselves we plotted the 
differences between the predictions at truncation level fc and fc + 2 on a log-scale. Fig. [3] (a) shows 
the change of (crjcrf). Apart from a cycle with period 2 (which is a peculiarity of the Neel state, 
see below) we clearly see the exponential decay of the differences. Fig. [3J (b) shows the changes of 
(erf erf ). Here we see the same behavior with approximately the same exponent, except for the points 
at k = 2. The fact that the truncated GGE predictions at k — 2 do not fit the exponential decay 
is not surprising: it is expected that for local operators spanning j sites the first j charges need to 
be added to reach the GGE prediction or at least the region where exponential convergence sets in. 
Similar behaviour was observed in the Ising chain in |13| . 

• At A = 2 and A = 3 the numerical data does not show convergence in the regime k < 12. In fact, 
plotting the differences on a log scale as before results in seemingly random data (Figs. [3] (c) and 
(d)). The correlations are monotonic functions of k and they are bounded so eventually they have 
to converge. We conjecture that there is exponential convergence also for A = 2 and A = 3, but it 
starts at some k — k c with k c > 12. 

• For large A all GGE predictions are close to the Neel values, where the z — z correlations are maximal 
and the x — x correlations are zero. For smaller A we see gradual departure towards the isotropic 
limit. 

Summarizing the above findings the following picture emerges. There exists a "convergence length" k 
which depends on A, such that for every short-range correlator O 



(0) (k) = (0) GGE +(a + P{-l) k ' 2 ) e- k '* + . . 



for k > k r 



(5.2) 



where a, f3 and k c depend of 0, and if O spans at most j sites then k c > j . The convergence length is 
small when the pre-quench and post-quench Hamiltonian are "close", in our case when A is large. The 
dots on the r.h.s. of (I5.2[) denote contributions which decay faster than the leading correction. The form 
of (15. 2p takes into account the observed oscillation with period 2, which is a peculiarity of the Neel state. 

Due to the limitations of our data we do not perform any k — > 00 fits of the correlators. For any 
practical purposes we suggest to take the values at k = 12. For larger A these data are already quite 
accurate, whereas for smaller A they are to be understood as lower or upper bounds, depending on the 
correlator in question. All our data shows clear difference between the thermal prediction (fc = 2) and the 
approximation of the full GGE (fc = 12). The differences are more pronounced for the 3-site and 4-site 
correlators. 

We conclude this section by providing a possible explanation for the observed cycle of period 2 in Fig [3] 
(a) and (b). For any even fc the charge Qk has terms proportional to (cr|<7L_ m — 1) with m = 1 . . . fc/2 |31j . 
where we already subtracted the ferromagnetic expectation value. In the Neel state these terms evaluate 
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to (— l) m — 1; they vanish for every even m. Therefore, whenever k = Am, the last two-site z — z operator 
added to the GGE does not feel the antiferromagnetic order of the Neel state, in contrast to the charges 
k = Am + 2 which are the ones who "submit" this information. This also explains the peculiar behavior of 

(k) 

the thermal multiplier p\ : its value changes considerably when a new charge with k = Am + 2 is added, 
but it almost stays the same when the truncation level is raised to k = Am + 4. 

6 Conclusions 

We considered global quantum quenches in the XXZ spin chain and showed that the Generalized Gibbs 
Ensemble can be implemented within the Quantum Transfer Matrix framework, and it yields predictions 
for the long-time limit of short-range correlators. The central idea was to construct truncated GGE's with 
only a finite number of higher charges; we then investigated the convergence of the predictions as the 
truncation level is increased. 

In our example we considered quench dynamics starting from the Neel state. Quenching to a finite but 
large A we found that our iterative procedure converges exponentially fast as a function of the truncation 
level. On the other hand, for smaller A convergence was not yet reached in the truncated GGE's with the 
first 12 charges. 

Correlations in the long-time limit are close to their values in the Neel state when we quench to a 
large A. On the other hand, quenching to smaller A the system departs towards the isotropic correlations. 
An important result of the present work is that the addition of the higher charges forces the correlations 
towards their Neel values. This is consistent with the general idea behind GGE, namely that the existence 
of the higher charges constrains the dynamics and physical observables can not relax to their thermal 
values. 

One of our original expectations was to find that adding k charges to the GGE would fix all correlators 
spanning at most k sites, possibly with negligible correction terms. Our findings show that this is only 
true, when the pre-quench and post-quench Hamiltonians are sufficiently close. Furthermore, even in these 
cases the dependence on the truncation level looks very similar for all correlators (see fig. H]). In our case 
we found that the x-x correlations reveal the expected pattern: exponential decay of the corrections sets 
in when the required number of higher charges is already added (compare figs.^a) and^b)). 

Our results show that the convergence length k (governing the convergence as we add the higher 
charges) is small for "small quenches" and that it is a new length scale which is independent of the 
correlation lenghts of the system. It is an interesting question how k depends on the pre-quench and 
post-quench Hamiltonians in more general situations. 

It would be interesting to consider the quench from the Neel state at the isotropic point (A = 1). 
Although this particular case is expected to exhibit the slowest convergence as a function of the truncation 
level, one advantage would be that in the XXX case all terms in the higher charges are known analytically 
|30) . therefore it might be possible to reach much higher truncation levels. 

The present methods could be applied to other quench situations as well. The only requirement is that 
the mean values of the conserved charges in the initial state should be evaluated exactly. In principle this is 
possible even for a finite interaction quench A — ¥ A, because the charges of the post-quench Hamiltonian 
are nothing more than certain combinations of Pauli matrices, and their mean values in the pre-quench 
ground state could be obtained by already available techniques. 

We would like to stress that in this work we did not attempt to prove the GGE hypothesis. Instead, 
our numerical results could be used as test of the GGE. We determined both the thermal and GGE 
predictions, and we expect that independent numerical investigations could confirm one of the two, or 
they could point to a different result. Numerical methods which simulate real-time dynamics typically 
involve a certain amount of integrability breaking, therefore they are expected to converge to our thermal 
predictions (the rows with k = 2 in the Tables 1 and 2), possibly with a pre-thermalization regime where 
the GGE applies. 
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Table 1: The Lagrange multipliers of the truncated GGE's and the predictions for the short-range correlators 
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Table 2: The Lagrange multipliers of the truncated GGE's and the predictions for the short-range correlators 
as a function of the truncation level 



19 




10 12 14 





■ 












0.2 




■ 


■ 


■ 


■ 


■ 


0.15 












- 


0.1 












- 


0.05 












- 




A 


A 


A 


A 


A 


A 




• 


t 


• 


• 


— •— 


— • i 



2 4 6 



10 12 14 



(a) |<7i 02,3,4 1 at A = 5 



(b) |afaf i3 ,4|at A = 5 




10 12 14 





■ 












0.25 




■ 








- 








■ 


■ 


■ 


■ 


0.2 












" 


0.15 












- 


0.1 












- 


0.05 
















A 


A 


A 


A 


A 


A 




• 


• 


• 


• 


• 


• 



2 4 6 



10 12 14 



(d) K crf,3,4| at A = 4 





■ 














0.3 




■ 


















■ 


■ 


■ 


■ 




0.25 














" 


0.2 














- 


0.15 














" 


0.1 














- 


0.05 


A 
• 


A 
• 


A 


A 


A 


A 


- 




, 


, 


• 


t 


t 


• 





2 4 6 



10 12 14 



(c) \a\ 02,3,4 1 at A = 3 



(f) K<rf,3, 4 |at A = 3 



Figure 1: The short-range correlators predicted by the truncated GGE's as a function of the truncation level. 
The top, middle and bottom curves correspond to o"i<t" +j - with j = 1,2, and 3, respectively. We plotted the 
magnitude of the correlators, the signs are given by (— 1) J . 



20 









■ 


i 


■ 






0.6 


■ 


■ 










- 


0.55 


- 












- 


0.5 


- 












- 


0.45 


- 












- 


0.4 


- 




A 


A 


A 


A 


- 


0.35 


A 


A 












0.3 


- 












- 


0.25 


• 


• 


• 


• 


• 


• 


- 



2 4 6 



10 12 14 





■ 


■ 












0.35 






■ 


■ 


■ 


■ 


- 


0.3 














- 


0.25 














" 


0.2 














- 


0.15 














- 


0.1 


A 


A 


A 


A 


A 


A 


" 


0.05 


• 


• 


• 


• 


• 


• 





2 4 6 



10 12 14 



(a) |crfcrf 3 4 | at A 



(b) kf(7f,3,4| at A = 2 



Figure 2: The short-range correlators predicted by the truncated GGE's as a function of the truncation level. 
The top, middle and bottom curves correspond to o"i<t" +j - with j = 1,2, and 3, respectively. We plotted the 
magnitude of the correlators, the signs are given by (— 1) J . 
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Figure 3: log(|0( fe+2 ) - 0( fc )|) as a function of k. 
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